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Abstract. We consider a region M in K" with boundary dM and a metric g on M conformal 
to the Euclidean metric. We analyze the inverse problem, originally formulated by Dix [5], of recon- 
structing g from boundary measurements associated with the single scattering of seismic waves in 
this region. In our formulation the measurements determine the shape operator of wavefronts outside 
of M originating at diffraction points within M . We develop an explicit reconstruction procedure 
which consists of two steps. In the first step we reconstruct the directional curvatures and the metric 
in what are essentially Riemmanian normal coordinates; in the second step we develop a conversion 
to Cartesian coordinates. We admit the presence of conjugate points. In dimension n > 3 both steps 
involve the solution of a system of ordinary differential equations. In dimension n = 2 the same is 
true for the first step, but the second step requires the solution of a Cauchy problem for an elliptic 
operator which is unstable in general. The first step of the procedure applies for general metrics. 

1. Introduction. We consider a region, M, in M. n with a smooth boundary 
dM. We assume that there is a Riemannian metric, g, on M that is conformal to the 
Euclidean metric with conformal factor v~ 2 where v £ C°°(M) is strictly positive. 
This means that g(x) = v~ 2 e, where e is the Euclidean metric, or, in Cartesian 
coordinates x = (x 1 , . . . ,x n ), g iJ (x) = v(x) 2 5' l: > . We analyze the inverse problem of 
reconstructing g based on measurements of the curvature of wavefronts produced by 
point diffractors located inside M and propagated according to the wave operator 
on (M, g). Indeed, geodesies for the metric g are rays following the propagation of 
singularities by a parametrix corresponding to this wave operator. In the seismic 
context v is the wave speed. As originally formulated by Dix [9], this type of data 
may in some cases be reconstructed from reflection data by variation of source and 
receiver locations at the surface of the earth. In particular, it is possible to recover 
from reflection data the shape operator for the wave front produced by a given point 
diffractor where the rays beginning at the diffractor intersect dM orthogonally [15 . 
Dix developed a procedure, with a formula, for reconstructing one-dimensional wave 
speed profiles in a half space from reflection data. Since Dix various adaptations have 
been considered to admit more general wave speed functions in a half space. We 
mention the work of Shah |23| . Hubral & Krey [13], Dubose, Jr. [TO], and Mann [T5] . 
We consider here the case of higher dimensional regions with Riemannian metrics 
conformal to the Euclidean metric. Our method is different in the cases of n = 2 
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and n > 3 and in fact we expect better results in the case n > 3. The problem is 
closely related to the problem where broken geodesies are observed on the boundary 
[IB] or when the Cauchy data of the solution of the wave equations are observed on 
the boundary, see e.g. [31 HI [T5] and references therein. 

Assuming that we know v and all of its derivatives on dM (c.f. [IE]), we may 
extend v to a function, which we will also denote by v, on a complete manifold M 
compactly containing M (in fact we can take M = M. n ). The corresponding extended 
metric is also denoted simply be g. As described in detail below, we measure the 
curvature of generalized^spheres for g centered at "diffraction" points intersected 
with an open subset of M \ M. From these data, and assuming we know v in M \ M, 
we show an explicit method to determine the function v in the Cartesian coordinates 
x = (x , . . . ,x n ) along geodesies of g which connect the diffractions points to the 
measurement region. This method can be viewed as a generalization of the work of 
Iversen and Tygel [Hj. We now proceed to introduce the concepts and notations 
necessary for the statement of our main results. 

For any (x,r)) £ VIM (fi indicates the unit sphere bundle with respect to g) we 
will write j x ^ v for the geodesic with initial data 7 a ,^(0) = x, j x ,n(ty — V- F° r r m the 
domain of 7^ we let C r (x, 77) denote the set of times i such that t is conjugate to r 
along "f XtV (by this we mean that (t — r)^ xri (r) is a critical point for exp 7 ^ ^ r \). 

For the moment we fix (xo, rjo) £ Cl(M\M) and use the notation C r for C r (xo, tjq). 
We will refer to the image of the set 

U G T y M : |£| s = R} 

under the exponential map as the generalized sphere of radius R centered at y. When 
t > r > is in the domain of "/x a ,r) an d t (£ C r there is a small portion of the 
generalized sphere of radius t — r centered at yt := Jxo,vo(^) containing J Xo ,vo( r ) that 
is an embedded submanifold S r t of M, Indeed, in this case we can define a vector field 
y r< t in a neighborhood of Jxo^oi 7 ") by writing for £ £ T yt M in a small neighborhood 
of -(t-r)% , m {t): 

s=l 

Geometrically, v r t gives the outward pointing normal vector fields to a part of the gen- 
eralized sphere centered at yt near r y Xo ,ria( r )- The shape operator S r .t £ (T 1 1 ) 7xo 
of T, r j at 7x ,j7 (r) is then given by 

for all X £ T 7xo ^^M, where V is the Levi-Civita connection for g. For the recon- 
struction of v we assume that So.t is known for alH > such that t ^ Co. In reflection 
seismology one refers to S rjt as the (partial) front of a point diffractor located at yt- 
We now introduce a mapping which defines local coordinates in which we will 
perform our initial calculations. We begin with picking a large to > in the domain 
of 7x , no such that to ^ Co. Next, let us take local coordinates x = (x 1 , . . . ,x n " 1 ) on 
S ,t snch that (x 1 , . . . ,x" _1 ) = defines x and suppose $ to : S ,t ^ U C 
is the coordinate map. We will assume that the image of these maps, U, is the same 
for all to- For x £ U, let 7!° be the geodesic with a special choice of initial data: 



^r,t(exp yt (C)) 



ds 
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Fig. 1.1. An illustration of the various notations introduced. 



7jr°(0) = — h>Q tQ {<& t 1 (x)). Then we define coordinates, (x,r), on some set by the 
inverse of the map 

*t (%r)=4°(r). (1.1) 

We define 

W(to) = *t {{(x,r) eUxR : r < t , r £ C to (^(2), -Vo tto ($£@)))}) 

and so that the map *S>t defines a local parametrization on W(to). It may not be 
a global parametrization because it may not be injective. The (x, r) local coordi- 
nates on W(t ), basically, are Riemannian normal coordinates centered at y to , but 
parametrized in a particular wayj^J? can be thought of as a parametrization of part 
of the sphere of radius to in T Vt M , and then r corresponds to the radial variable in 



T yto M. Figure 1.2 shows a possible depiction of W(t ) and the coordinates defined 
there. We note that the domain W(t ) includes "fx ,ri ([0, to)) \ {lx ,r l0 ( r ) '■ r € C to }. 
We also define for some L > 

W= (J W(t ). (1.2) 
t e(o,i) 

In our reconstruction we cover M by sets like W(to), and then recover v on each 
of these regions separately. In the case that there are conjugate points to t along 
1x0,7)0 we wm a l so have to cover some regions with more than one such set. Note that 
along 7x ,?7o t ne coordinate vectors d/dx^ are Jacobi fields, and are defined even at 
the conjugate points. 

Finally, we introduce frames {Fj°(x,r)}% =1 defined by parallel translation along 
7g° such that for j — 1, ... , n 



F°(x,0) 



(£,o) 



where we are using the notation r — x n and so 

F*o($,0)= 7 | , (0) 
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Fig. 1.2. A figure illustrating the coordinates (x,r) defined on the set W(to). 



points in the opposite direction of vot {^t (x, 0)). To simplify the presentation we 
adopt this notation, r = x n , throughout the paper. Also we write {/j? (x, r)}J =1 for 
the corresponding dual frame; that is 

(fi(%r),F t k °(x,r))=6l. 

where (., .) denotes the usual pairing of T t 0( ,M and T* t0 M. In the sequel will 

also consider the shape operators <SV,t when xq is replaced in the above construction 
by another point in Eo.t represented in the coordinates by x. We thus have for 
each x and < r < t < t such that r and t are not conjugate along "y~.° a tensor 

Sl° t (x) € (Ti) tp,-.M. We represent Sl° t (x) using the frames constructed above as 

Sj? t (x) = st(x, t ; r,t)fi (x,r) ® ^"(^r). (1.3) 

For fixed to and x we will also use the notation Sj(r, t) — Sj(0,to;r,t). Note that 
immediately from the definition we have s"(x, to;r,t) = s 3 n (x,t ;r,t) = for all j 
and because of this in what follows when we write s(x, to ;r,t) (respectively s(r, t)) 
without indices we will actually be referring to the (n— 1) X (n— 1) matrix s^(x, to; r, t) 
(respectively s^(r, t)) with j, k = 1, ... ,n— 1. The data for our recovery are the matrix 
elements s*?(2:, to; 0> t), and their first three derivatives with respect to t, for < t < to 
and t not conjugate to along 

In [8], we obtain the following result: It is possible to uniquely determine the 
Riemannian metric g in a neighborhood of 7x ,n D ([0,to)) m Ricmannian normal co- 
ordinates having origin at the point yt a (this can be done for general metrics, not 
just ones which are conformally Euclidean). Here, we cast this result into an al- 
gorithm, and construct a conversion from the mentioned coordinates to Cartesian 
coordinates, which is the main contribution of this paper. Essentially, we generalize 
the time-to-depth conversion in Dix' original method to multi-dimensional manifolds 
with Riemannian metrics conformal to the Euclidean metric and, roughly speaking, 
show that if we measure near the point xo the shape operators of the wave fronts of 
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waves diffracted from the points 7x0,7)0(^0), we can then determine the wave speed v 
near the geodesic "fx ,vo- The theoretical contributions of this work are thus contained 
in the following two theorems. 

Theorem 1.1. Suppose that M C M — K n and M \ M is op en and nonempty, 



v 6 C°°{M), W is a set of the form constructed above (see [1-2)) using the metric 
g = v{x)~ 2 5ijdx l Ax^ , and So,t C M \ M for all to G (0,L). We also assume that 



(M,g) is complete. If v \^M\M)nw * s known then from the data s(x,to;0,t) (see (1.3)) 
for all tg in an open subinterval of I — (Ti,T2) C (0, L) n Co(xo, rjo) c , x G U, and 
t G (0, to) H Ct ($^ ) 1 (S;), — Vo,t {^t 1 (%))) c it * s possible to recover v in a neighborhood 
of 7rc ,j) ([0? ^2)) in Cartesian coordinates. In dimension strictly larger than 2 the 
reconstruction only involves solving ordinary differential equations along the rays of 
9- 

Note that the sets C r (x,n) are always discrete, and so the hypotheses of the 
theorem assume we have data except for at a discrete set of times. We also stress 
here that the reconstruction is local along geodesies. That is, to reconstruct v in a 
neighborhood of 7^0,770 we only require measurements of the shape operator near the 
point xq for generalized spheres £o,t centered at points 7z , I?0 (io) with radii to > 0. 
We recall s(x,to;Q,t) is representation of the shape operator of these generalized 
spheres in the local coordinates and that we have xq G £0^0 • ^ ne case that the 
dimension n > 3 the conversion from Riemannian normal to Cartesian coordinates 
involves solving a system of n + 3n 2 + n 3 nonlinear ordinary differential equations. In 
the two dimensional case the construction requires the solution of a Cauchy problem 



for an elliptic operator (i.e. we must solve the scalar curvature equation (2.5| with 
known boundary data). The discretization of the system is directly related to the 
available "density" of scatterers. 

In orderjto state our second result, we musl^ntroduce generalized distance func- 
tions. As M is complete, the map exp y : T y M — > M is surjective, and by Sard's 

theorem, the set C(y) C M of critical values of exp y has measure zero. Suppose 
x G M\C(y) and let £ G T y M be such that exp y (£) = x. Then £ has a neighborhood 
Vy^ G T y M such that exp y : V Vt £ — > V y ^ = exp y (V y ^) is a diffeomorphism; one says 
that exp" 1 : V' * — > V y ^ is a local inverse of exp y corresponding to £. In V y ^, we 
define the function, p(-;y, £), by 

p{z;y,0 = lexp-^z)^, 

and call this function the generalized distance or travel time function from the point 
y associated to direction £ ; the wave fronts observed from a point source at y give us 
its level sets. We may now state our second theorem. 



1.1 



and suppose that dM is 



Theorem 1.2. J^et M and M be as in Theorem 
smooth and W C M is an open set such that for all y G W f~l M i 
non-normalized geodesic 7 for g such that 7([0, 1]) C W , 7(0) = y, 7(1) G dM and 
7(1) (fz. TdM . We use the notation V = W' D dM , and also assume that w li^'n(M\ m) 
is known. Further suppose that A is an indexing set such that there is a bijective map 
from A to 

{(y,0€T(W'nM) : V^nT^V}, 



() 



and label y, £) = p\ according to this map. Then from knowledge of p\\r for A 6 A 
we can recover the wave speed v in W' . 

This generalizes an earlier result which says that the boundary distance func- 
tions of a compact Riemmanian manifold given on the whole boundary determine the 
manifold uniquely, see [T7] and [THl Section 3.8]. 

The structure of the rest of the paper is as follows. In section [2] we go over 
some background from Riemannian geometry that is necessary. Section [3] reviews 
the first step of the recovery procedure in which we reconstruct the metric g in the 
coordinates (x,r) given by the map ^ to . Then sections [4] and [5] describe the second 
step of the recovery procedure in which the geodesies and wave speed are reconstructed 
in Cartesian coordinates respectively in the case of three and higher dimensions and 
the case of two dimensions. Finally, the proofs of Theorems |1.1| and |1.2| are given in 
section [6] and section [7] contains some concluding remarks. We also provide detail on 
the conversion from travel time measurements to measurements of the shape operator 
of generalized spheres for the case when dM is flat and v is constant in a neighborhood 
of dM in E 

2. Preliminaries. We summarize the basic differential equations from Rieman- 
nian geometry that we will use. We mention some general references to Riemannian 
geometry [TTJ [2T] [52]. In this work we will use the conventions from [2T] for the 
curvature tensor and related quantities in local coordinates. 

2.1. Geodesies. We evaluate the geodesies by solving 

dV ■n,(7(*))^ = I (2.1) 



At 2 wwv " At df 

where 



r M - 1 nP i f d9kp + d9lp d9kl 



are the Christoffel symbols. It is also possible to find the solutions of (2.1 ) using the 




Hamiltonian flow for the Hamiltonian H(x,p) = \pjp^g : ' (x). Although we will not 
use the Hamiltonian formulation here, we note that it gives the system 

- -A p .p Qg^W) 

for the geodesies. From this, we see that, in terms of seismic ray tracing, the geodesies 
may be identified with generalized image rays. 

In our case, assuming isotropy (i.e. that the metric is conformally Euclidean), we 
have 

Km = Wn + - S qm 6 M ) / = l0g( V ). 

It is more convenient to work with / than v and we will generally do so throughout 
the remainder of the paper although it is clear that recovery of / is equivalent to 
recovery of v. To make the notation more concise below we introduce the shorthand 

^qrn ~ u q u rn t <J m <J q u qm u . 
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2.2. A frame, parallel transport. As mentioned above, Fj°(x,r) denotes the 
parallel translation of djdx? along 7!° from to r. Thus, for every x and r > 0, 
{.F-f (x, r) , . . . , (x,r)} forms a basis for T^ M. The invariant formula for parallel 
translation is 

V^ (r) F*°(z,r) = 0. 

If we introduce matrices which give the frames F?j (x, r) in Cartesian coordinates 

d 



then the invariant formula implies that these satisfy 

^ r Fl + {il°) k {r)T l km F? = Q 



(2.2) 



-FUt ;x,r) + F*(t ;x,r)rlM S (r))F™(t ;x,r) = 0, (2.3) 



or 

dr 3 ' 

and since the coordinates x on So,t are known, we also know the initial conditions, 
Fj(* o ;£,0). 

2.3. Curvature. The Riemannian curvature tensor in any coordinate system is 
given as a (1,3) tensor field by 

nP _ _^_rP _ d r p r ,: r p _ r i r p 

3fe/ ~~ feJ fc/ dx k i l kl i { jl ki ' 



or as a (0, 4) tensor field as 



jkl — Rjklp- 



The Ricci curvature tensor is given by the trace 

Ricij = R^j 

and the scalar curvature is seal = g lJ RiCij. In the Cartesian coordinates the Riem- 
manian and Ricci curvature tensors are given respectively by the following formulae 
in terms of /: 

R = e- 2 f (5 (-Hess(/) - V/ ■ V/ + ||V/| 2 ^ 

where is the Kulkarni-Nomizu product (in coordinates, see [21j for the invariant 
formula), Hess(/) is the Hessian matrix of second derivatives, V/ is the (Euclidean) 
gradient, and |V/| is the Euclidean norm of V/; 

Ric = (n - 2) (Hess(/) + V/ ■ V/) + 5 V] (a/ + (2 - n) |V/| 2 ) . (2.4) 

Also, when the dimension is n = 2 the scalar curvature satisfies the so-called scalar 
curvature equation 

A g f = l - seal (2.5) 
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where A g is the Laplace-Beltrami operator corresponding to g given in any coordinate 
system {y j }]=i by 

Here g = det([gjk] 2 _ k=1 ) and [g^] 2 f-=i IS the inverse matrix of [gjk] 2 : k=i- From these 
formulae we see a fundamental difference between the two dimensional case and the 
case of three or more dimensions. In two dimensions the Ricci curvature and scalar 
curvature give only the Laplacian of /, and so to find / from these curvature tensors 
would require the solution of an elliptic equation. On the other hand, in three or more 
dimensions the Ricci curvature tensor depends on all the second partial derivatives of 
/, and in general we can recover a formula for Hess(/) in terms of the Ricci curvature 
and V/. Indeed, if we define 

G = ffic~(n-2)(v/-V/-|V/| 2 <5) (2.6) 



then from (2.4) we may calculate 



Hess(/) = -l^( G -^-L_tr (G )^. (2.7) 

This is possible in three dimensions, but not in two dimensions, and is the reason we 
must consider the two cases separately. 

We will also write the Riemannian curvature on the geodesic 7-° in the frame 
obtained by parallel transport as 

i?yfo;z,r)F*°(x,r) = i? 7 , n(r) (F;°(^r),F*°(S;,r))^(3?,r), 



R%l (to \x,r) = (fi (x, r) , (r) (if (x, r) , F*> (x, r ) )Fp (x,r)). 
Recalling that F*°(x,r) = 7- (r), we also write 

^ P j(to;x,r) = R^ nn (t ;x,r), 

and for fixed t 

r?(r) :=r^ o ;0,r) = (/£ (0, r), R^^F^O, r), % „ {r))% 0tVo (r)) , (2.8) 

for the directional curvature operator which we reconstruct in the first step of our 
procedure. Note that as with s, for any j r™(r) = r J n (r) — 0, and so when we write 
r without indices we will actually be referring to the corresponding (n — 1) x (n — 1) 
matrix. 

We continue in the next sections to describe the actual reconstruction algorithm. 

3. Reconstruction procedure Step 1: Determination of the metric 
in (x, r) coordinates. The reconstruction procedure consists of two steps. In the 
first step we consider only the single geodesic ^xa^o an( i reconstruct r^(r), and then 
the metric g as a function of r for x = 0. In the second step, we determine /, and 
therefore also u, in a neighborhood of 7a; ,?j by also varying x and <o- 

Following the geometric analysis of [S] we now describe the first step of the pro- 
cedure which is itself broken up into a number of substeps below. 
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1. Let W = V 3 (r, t), j = 0, ... ,3 represent (n — 1) x [n — 1) matrices. We 
solve the autonomous system of ordinary differential equations for V(r,t) — 
(V J >,*))f=o, 

/- ±V°(TV 3 )V°, 
^(V^TV^V + V"(TV 3 )V 1 ), 

±(V 2 (TV 3 )V° + V°(TV 3 )V 2 + 2V 1 (TV 3 )V 1 ), (^-l) 
i(V 3 (TV 3 )V° + V°(TV 3 )V 3 + 3V 2 (TV 3 )V 1 + ^{TV 3 ^ 2 ), 



£v° = - 

i-y 3 = - 



in which 



(TV 3 )(r) = V 3 (r,r), 



for < r < t < to. This system is supplemented with initial data, 

W(o,t) = {af(s(o,t))- 1 } J 3 =0 . 



(3.2) 



Applying Picard's theorem in a standard way, we see that the equations (12)- 
(13) have unique solution on some interval t € [0, t{\ (for details on this 
see [8]). In practice, we may use a Runge-Kutta method to solve the system 
numerically for < r < t < t\. The system will not generally have a solution 
all the way up to to', in this case, we must divide the interval [0, to] into several 
subintervals, and reconstruct on each of these in turn as described below in 
substep 3. 

2. We extract the directional curvature operator, 



r(r) = \(TV)(r) 



(3.3) 



Note that this matrix r(r) is (n — 1) x (n — 1), but recall that from this we 
can recover the full directional curvature operator r^(r) since the nth row 
and column of r^(r) are both equal to zero. 
3. In general the first two steps only reconstruct r^(r) for < r < ti where 
t\ < to since we may not able to solve (3.1) all the way up to to- Here, we 
describe how to find rl(r) for r in the entire interval [0, to]- The idea is to 
replace by t\, and then use knowledge of Sj(t\,t) for t>t\. Indeed, let us 
fix t > t\. Now we find the coordinates of the Jacobi fields corresponding to 
the coordinate vectors fields of (x, r) along 7 XOir;o with respect to the parallel 
frame. These can be found by solving the system 



d_ 
Or 



JIM) 



6$ 
-rj(r) 



(3.4) 



with initial conditions 



mt) 



si 



(3.5) 



We can then recover s(ti,t) by the equation 
dj 



Or 



{h,t)(r 1 )(ti,t) = -s{t 1 ,t). 
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Region 1 



Line 2 



Line 1 



Fig. 3.1. 



A depiction of how the first step of the algorithm pr oceed s. The original data give s 
on line 1. Then substep 1 recovers {V^}|__ j in region 1 by solving 3.1). Next substep 2 gives r(r) 
for < r < t\ by K3. 31), and then substep 3 gives s on line 2. Returning to substep 1 again gives 



{V J }| =1 in region 2. Then as before substep 2 gives r(r) for t\ < r < t2 t and substep 3 gives s on 
line 3. Continuing in this way we reconstruct r as far along 7^ ,rj as we like. 



Now we may return to substep 1 and solve (3.1) with (3.2 1 replaced by 



Vi(t 1 ,t) = {6i( S (t 1 ,t))- 1 }* 



This then allows us to use ( |3.3[ ) to recover r(r) for r up another time t 2 say 
with ti < t 2 < to- If ti is less than t we repeat this same procedure again 
with t\ replaced by t 2 , and so on. By results in [5] there is a lower bound 
on the size of each step we take (i.e. a lower bound on U — and so by 

induction we eventually recover r^.(r) for r on the entire interval [0, to]. For 
a visual depiction of how the first three substeps proceed see Figure [3] 
We now obtain the metric in the (x, r) coordinates given by the coordinate 
map ^ to along 7a; 7i , which we write as p_yfe(0, r), by the formula 

9jk(0,r) = j^(r,t )j^(r,t )g pg , (3.6) 

where 

g pq = g(F p (0,0),F q (0,Q)) = v(x )Fi(0,Q)F l q (0,0)5. l3 . 



We assume in the hypotheses of Theorem 1.1 that v(xq) is known, and F^ (0, 0) 
(resp. Fjj(0, 0)) are the components of the coordinate vector d/dx p (resp. 
d/dx q ) with respect to Cartesian coordinates. Thus g pq is known. By ad- 
justing the choice of xq we also find the metric with respect to the (x, r) 
coordinates where they are defined (i.e. on all of W(to)). 

4. Reconstruction procedure Step 2: Transformation of coordinates. 

By the procedure described in the previous section we can reconstruct the metric 
gjk{x,r) in the coordinates (x,r) everywhere in the domain W(t ) where those coor- 
dinates are defined. Thus we can also reconstruct the Ricci curvature tensor in these 
coordinates and also the scalar curvature, which we will use below. In this section, 
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we show how to determine the velocity function w(7« a (r)) and the geodesies 7~°(f) in 
the Cartesian coordinates from this information. The reconstruction is done initially 
up to the first conjugate point to to along 7x ,J7 - Then, as described at the end of 
this section, to must be changed to allow reconstruction past this conjugate point. 

As observed above, the coordinate vectors in the (x, r) coordinates are Jacobi 
fields along 7!°; in particular, if 



d 

dx k 



= j l k (t ; x, r)Ff° (x, r) = j{(t ; x, r)F?(t ; x, r) 



7?W 



d 
dxP 



(recall that are the coordinate vectors for the Euclidean coordinates) then the 
matrix ] l k (to; x, r) satisfies 



d_ f $ k {t ;x,r) 
dr \ j l k (t ;x,r) 

supplemented with the initial data 



j l k (t ;x,0) 
i[(to;x,0) 



8' p 
-r l (t ;x,r) 



(4.1) 



61 



-s l k (x, t o ;0,t o ) 



In fact, here we are simply adding the dependence on x and to to the same quantities 
already considered in the previous section. Since parallel translation preserves the 
metric, we also have the relation 

v(l~° (0)r 2 Fj(t ; x, 0)F«(to; x, 0)S lq = «( 7 ^(r))- 2 Fj(i ; x, r)F q k (t ; x, r)S lq . 

By taking determinants of the matrices on each side and then the natural log, we 
obtain the following formula 



/( 7 |°(r)) = ilog(«(7| , (0)) n 



det(F l At ;x,r)) 



det(Fj(t o ;x,0)) 



(4.2) 



Recall that Fj(to;x, r) are the components of the parallel frame with respect to the 
Cartesian coordinate vectors (cf. (2.2)). Also, combining some of the previous formu- 
las we have 



^/(7| (r))=^(to^,r-)^(t ;x,r)g(7|°(r)). 
Away from conjugate points along 7!°, we may invert to obtain 

g( 7 |°W) = (F- 1 )f(io;S,r)(j- 1 ^(io;^r)^/(74 (r)) 



(4.3) 



(4.4) 



which may further be combined with (4.2 1 to show 



^(7l°W) = i(f- 1 )f(* ;$,r)a- 1 )i(to;£,r) 



r)rr to r)F b 

^-(x) + (F^)UUy,x,r)^(t ^,r) 



(4.5) 



where we use the notation 

a to (x) := log 



vhl°(o)Y 



|det(Fj(t ;x, 0))| 
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Remark 4.1. A "stepping" recovery procedure in the spirit of Dix' original 
method may now be presented. We introduce a step size h in r, and for a £ N we 
label r a = ah. If we know f(ji° (r _i)), Fj(t ;s, r a _i), and j^(t ;a:,r a _i), then we 
approximate the same quantities at r Q by the following strategy. First we numerically 
estimate the derivatives 



Then we use (4.4 1 to estimate the derivatives 



0). 



We then have an estima te o f F l km (-/^ (r Q _i)). Next, we use this estimate to perform 
a forward Euler step in (2.3) and get an approximation of Fj(t ; x, r a ). Then, finally, 



we use (4.2) to obtain the approximation for f{"i~{r a )). We note that 3 k (to;x,r a ) 
may be obtained through a completely independent calculation using (4.1). Because 
F^(t ;x,r) — ( A /~ ) l (r), we can then approximate (7- ) l {r a ) as well. 

A closed system of ordinary differential equations for n > 3. While the 
technique described in remark [4. 1| may be a direct generalization of Dix' method, in 
fact we seek a single closed system of ordinary differential equations that could be 
solved using any numerical scheme to give all the desired quantities. This is possible, 
although the method we present here works in three or more dimensions only. The 



reason for this limitation, as explained above, comes from our use of formulas (2.6) 



and (2.7) to express the Hessian of / in terms of the Ricci curvature and the first 
order derivatives of / in Cartesian coordinates. Since we actually know the Ricci 
curvature in (x, r) coordinates given by ^t , which we will label as Ric pq , we also 
need the formula for the tensorial change 



(4.6) 



Now we can describe how to get the closed system of ordinary differential equations. 



We differentiate (2.3) and use (4.3) to get 



d_9F[ = dFl pmG ik df_ 
drdx a dx a 3 qm dx k 



dF] 



r)f F) 2 f 

n dx a qm dx k n 3 qmJa p dx k dx c ' 1 ' 



where we have suppressed the dependence on (t a ;x,r). We may use (2.6), (2.7), (4.5), 
and (4.6) to express the right-hand side of (4.7) only in terms of j, F, derivatives 



of F, and the Ricci curvature Ric pq . Combining all the previous equations we now 
have a closed system of ordinary differential equations that may be solved uniquely 
up to conjugate points. In the next paragraph, we summarize the entire method for 
convenience. 

As claimed above we have now produced a closed system of ordinary differential 
equations that may be solved to obtain v and 7^° in Cartesian coordinates. The 
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system is nonlinear and contains n + 3n 2 + n 3 equations. It may be written as 



d_ 

Or 











ii 




n 




^dxP ' 





dx ' 



Wj. k (r,x,j,j,F,% 



W l F . k (r,x,3,'j,F,%) 



(4.8) 



V 



We describe how each of the "W" functions on the right-hand side are to be evalu ated. 

and Wi. k are the simplest. Recalling that F^{x,r) = j^(r) and using (4.1) they 
are given by 



W\=F l n and Wl. k =] 



Next, again according to (|4.l|), W- is given by 



W! 



j;fc 



-r l \ j 



Since x n = r, W l F . k is given by 



wl. 



F:k 



dx 11 



Finally, W l aF , is given by (|4.7|) where we calculate 



§£j using (|4.5 



and calculate 



a ®iQ x j in several steps using the values of already calculated, d4.6|, (|2.6|, and 
then (2.7). System (4.8) gives a nonlinear system of ODEs which can be solved to 
recover 7 andi^ up to conjugate points. However, at the conjugate points the matrix 
j will become singular, and we will not be able to continue. 

Note that we do not explicitly solve for v in this s yste m, but after we have found 
F k , then / and therefore v may be calculated from (4.2). Indeed since / = log(u), 
(4.2) becomes 



det(FHt ;x,r)) 



det(Fj(t ;x,0)) 



l/n 



The presence of conjugate points. By its construction, the system (4.8) has 
a solution up to the first conjugate point of £o along "fx ,rj i an( i since the functions on 
the right hand side of (4.8) are Lipschitz continuous away from the conjugate points 
this solution is unique. Thus we can recover v along the entire length of 7 :Eo , ??0 using 



(4.8) if there is no conjugate point to to. However, if there is a conjugate point, then 
we must follow a more sophisticated strategy. 

Note that once we are able to calculate the matrix j(to;x,r) after step 1, then 
we can identify all of the conjugate points of to. Since we are only concerned with 
a finite length along the geodesic J Xo ,vo> ^ e Morse Index Theorem (see e.g. [20l 
Theorem 15.1]) there are only a finite number of points conjugate to to, and each of 
these conjugate points also has only a finite number of conjugate points on the finite 
interval of interest. Thus we can pick an alternate value of t' such that to and t 
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do not have any of the same conjugate points, and there are no conjugate points to 
either on [t , t ] . Performing step 1 with to replaced by by t' Q , we can also consider 
the system (4.8 1 with to replaced by t' . 

Now suppose that {tj}'jL 1 is the set of times conjugate to to in the interval [0, to) 
given in decreasing order (i.e. so that tj > tj+i for all j). Then for each j there is 
an open interval Ij C [0, t' ) containing tj such that Ij does not contain any times 
conjugate to t' along ^ Xo ,r lo - Suppose that {Ij}™=i is another collection of open 
intervals so that {ij}™^ 1 U {I'j}f =1 covers [0, to) and Ij C [tj, tj-x] for all j where we 
are setting t m+ i = —e for some e > sufficiently small. 

Now, let us set U m+ \ = U (recall that U is the domain of the inverse coordinate 



map $t ). By solving (4.8 1 and possibly shrinking U m+ i we can find the wave speed 
v in Cartesian coordinates on the set W m +i :— ^t {U m +i X Im+i)- By the continuity 
of ^ t ' w e can find an open set U' m C W 1 ^ 1 and an open interval J' m+ i containing 
and intersecting I' m nontrivially such ^t' {U' m x J' m+ i) C W m +i- Therefore, since 
the wave speed is known in W m +i we can find the parallel fields and Jacobi fields 
in Cartesian coordinates corresponding to t' Q for x £ U' m and r € J' m +\- Now, after 



possibly shrinking U' m , we can solve (4.8 ) with to replaced by t' to find the wave speed 
in Cartesian coordinates on ^f t > (U' m xl' m ). Thus we have reconstructed the wave speed 
on fy t > {U' m x (I m+ iU/^J). Switching the roles of to and t' we can reconstruct the wave 
speed on a set of the form ^ to (U m x (I m +i U I' m Ul m )) for some open set U m C U m+ \. 
Continuing in this way after a finite number of steps we can eventually reconstruct 
the wave speed on a set of the form fy to {Ui x (— e, to)) which is a neighborhood of the 
geodesic 7z 0j r )0 . This completes the reconstruction in dimension three or higher. 

5. Two-dimensional case. The method of the previous section, step 2 in the 
reconstruction procedure, will not work in two dimensions. The basic reason for this is 
that we cannot determine all of the second partial derivatives of / from the curvature 



of g, but rather can only obtain the Laplacian of / as shown in (2.5). Therefore in 



the two-dimensional case we are forced to recover / by solving (2.5). We discuss this 
in more detail below, but first we will also revisit step 1 in the two-dimensional case 
where the formulae can be simplified. 

The main simplification comes from the fact that in the two dimensional case 
the trace of Sl° t (x) contains all of the same information as >!?' t (x) itself, and so it is 
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actually easier to consider this trace. Indeed, let us define 

a(x,t ;r,t)=tr(S t ° t (x)). 



In this case we have the following simple method of calculating a(x, t ; r, s) away from 
conjugate points. If r is not conjugate to t along 7 ' 
defined for (z, s) in a neighborhood of (x, r) by 



~° , then there is a distance function 



exp 



??(*) 



In the seismic context this is nothing other than the local travel time along rays close 
to 7-° from 7~°(i) to points near J~(r). By 21, p. 46] we have 

a(x,t ;r,t) = A g d^' rt . 

We continue to review step 1 in the two dimensional case. 

5.1. Step 1 redux: The two dimensional case. We note that the {V J '}^_ 
which appear in equation (3.1 1 are actually all scalars and the only nonzero component 



of is which is what we recover in (3.3). Note also that in the two dimensional 



case, r} contains the same information as the sectional curvature (if F\ is chosen to 
have unit length with respect to g then in fact r\ is the sectional curvature). 

The other simplification occurs in substcp 4. In the two dimensional case the 
metric must have the form 



9jk = </°{x,r) 2 dx 2 + dr 2 
Now by [2TJ p. 46], ip to (x, r) satisfies the equation 



dr 2 



if to (x,r) +r{(t ,x,r)(p to (x,r) = 



(5.1) 



where v\(to, x, r) is already known. Since the metric is also known in a neighborhood 
of So, t we may simply solve this equation with the known initial data tp to (x, 0) and 
d r ip to (x,0) in order to find the metric in the (x, r) coordinates defined by ^ ta - We 
see that it is not necessary in this case to compute the matrix j corresponding to the 
Jacobi fields. 

Now we continue to show how step 2 may be accomplished in the two dimensional 
case. The method makes use of the scalar curvature rather than the Ricci curvature. 

5.2. Step 2 redux: The two dimensional case via the scalar curvature 
equation. In step 2 we must take a different strategy for the two dimensional case. 
The difference is that we cannot express the second derivatives of / which appear 
in (4.7) in terms of j, F, derivatives of F, and the Ricci curvature. Instead we use 



a method inspired by the treatment from |16[ Section 4.5.6] of a different problem. 
After we recover the metric g in the coordinates (x, r) given by ^t we use t ne scalar 



curvature equation ( 2.5 ) to directly solve for / in these coordinates. Indeed, we assume 
that 



v \zo, to and ||:|E , t0 



(5.2) 
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are known, and so in fact we have Cauchy data for / on So,t - Thus / satisfies a 



Cauchy problem for the elliptic operator A g (see ( 2.5 )) expressed in (x, r) coordinates. 
This equation is given explicitly by 

?(«d;5,r)- 1 / 2 ^?(to;x s r) 1 /V k (to;2,r)^/(T|'(r)) = \ scalar)). 



dx k 



Here gi (t ; x, r) is the inverse of the matrix given by formula (3.6 1 extended to values 
of x other than and g{to\x,r) is the determinant of the matrix gjkito) x, r). The 
scalar curvature, seal, can be computed from gj k (x,r). Expressed in the coordinates 
this is an elliptic equation although it is degenerate when r is conjugate to to along 
7~°. Adopting the notation of the previous section, by the unique continuation prin- 
cipal (for a modern review of Cauchy problems for elliptic operators see [2]) we can 
first reconstruct /(7~ (r)) for (x,r) G U m+ i x I m+ \. We note, however, that this 
reconstruction is generally unstable (once again see [5] for a detailed review of the 
stability of this type of problem) . 

Now once we have recovered / in (x, r) coordinates for (x, r) G U m +\ x I rn +i the 
system (4.8) can be replaced by a significantly simpler system. Indeed, if we combine 



the equation (4.1) for the Jacobi field matrix with (2.3) and (4.4), then we have a 
closed system of ordinary differential equations which may be solved just like (4.8 ) for 
the higher dimensional case. For convenience we write this system down explicitly. 
The system is 





(A 




d 






dr 













W{. k (r,x,i,'i,F) 
\W l F . k (r,x,j,'j,F)J 



(5.3) 



Here W l (r, x, j, j, F), WL(r,x,j,j,F), and 



formulae shown below (4.8) 

to ((231) and (F4|> b y 



(r,x,j,j, F) are given by the same 
The last entry on the right hand side is given, according 



W l F;k (r,x,3,i,F) = Q l * q F»F*(F- 1 ) i j (y 



dx 



M°{r)). 



Solving these equations we can recover /, and therefore also v, in Cartesian coor- 
dinates on ty to (Urn+i x 7 m _|_i). Once we have this we can find Cauchy data on the 
surface U' m x {t' m } for some t' m G I' m D {t < t m } for the scalar curvature equation in 
the coordinates given by ^ t ' , an d repeat the process described above. Thus we can 
introduce a stepping procedure as in the higher dimensional case and this completes 
the reconstruction in the case of two dimensions. 



6. Proofs. Proof. 1 1 . 1 
have presented in sections 



This theorem follows from the recovery procedure that we 
jj] through [5] Indeed, from results in [5] applying step 1, 
which is described in section [3j for any x G U we can recover the metric g in (x, r) 
coordinates corresponding to ^t f° r anv G !■ Then, in dimension 3 or higher, 
the argument of section [4] completes the proof, while in dimension 2 we must use the 
scalar curvature equation as described in section [5j □ 

The proof of Theorem |1.2| involves reduction to the case of Theorem |1.1| We also 
perform this reduction more explicitly in the case that dM is flat and v is constant 
in a neighborhood of dM in the appendix. 
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Proof. [L2] 

Let us begin by taking any A G A. By the hypotheses there exists a point 



z G domain(pA ) n T. Let x 



) be a set of local coordinates defined 



on an open subset of domain(pA ) ^ V containing z and denote by $ z the coordinate 
mapping which we suppose has range given by U C R n_1 . We construct $ z such that 
<& z (z) = 0. Finally, let v denote the outward pointing unit normal, with respect to g, 
vector field for dM. 

We will now write ddM for the exterior derivative on dM. For x G U and e > 
let us consider the set 



7v = {^A : daMpxi^^ix)) = d dM p\ ($ z 1 (x)), 
and Px (<5>- l {x)) < px^- 1 ^)) + e} 



(6.1) 



and mapping a £i£ : 72 



(-e,PA (^ 1 (^))) defined by 
AX)=p Xo (<p- 1 (x))-p x ^- 1 (x)). 



This set and mapping can be found from the data. Now, if range(a£ ;£ ) is dense in 
(— e,pA ($J 1 (a?)) for all x in a neighborhood of and \ddMP\ (^J 1 (x))\g < 1 then 
the geodesic represented by jo,e (see the next paragraph) is contained in W and 
intersects dM transversally. In this case we continue. Otherwise we do not use this 
choice of Aq. 

Intuitively, 75^ is a representation of the portion of a geodesic passing through 
$>~ 1 (x) which lies inside M, and the range of ctx,e parametrizes this geodesic segment. 
Indeed, for every A G 72. e , let V\(x) G T^-i^M be given by 



v\(x) 



a 



Px (sj 1 ^)) 



d 

dx k 



(6.2) 




dxi 



u^-\x)) 



9z j 



where g%. is the metric g restricted to dM expressed in the x coordinate frame. Then 
7x,e represents a segment of the geodesic 7^-1,^ _^ ,~. 

We now begin relating the constructions we have made so far with the data 
required to apply Theorem 1.2 First, let us take a sufficiently small constant e > 
and set 

£0 = 7z, VXo (o)(e), Vo = -7s> Ao (o)(e)' 

These can be constructed from the given data because Jz, Vx (o)((0,e)) C W C\(M\M) 
for e sufficiently small. We can also construct the map $ _1 : U — > M defined by 

Taking e sufficiently small, and possibly shrinking U we can make <i>~ 1 a diffeomor- 
phism onto its image, and so <fr is a coordinate map on the image of $ _1 which, using 
the notation introduced earlier in the paper, we take to be Sq t with i := Z+ p Xo (z). 
Also following the earlier notation we let vo,t 0?) be the outward pointing unit normal 
vector for So,t a f &~ 1 {x). 
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Suppose now that x £ U, and t £ [e + p\ Q {z) — p\ ($> z (%))t to]. If e + p\ a (z) — 
p^ ($J 1 (a?)) and t are not conjugate along 7$-i(gv_„ t ($) then there exists a A g 72 c 
such that p A (<I>J 1 (a;)) = i— e— p Ao (z)+p Ao ($J 1 (a;)). Now we have that Ua(-) is defined 
on a neighborhood [/' of ir, and the set S ,t is then given by 

Since we can recover this set and we know g £ M \ M we can calculate the shape 
operator S ° t (x) provided it exists. Therefore we may calculate s(x,t o ;0,t). Since 
there are only a finite number of t € [e + p\„(z) — p^^J 1 ^)), to] that are conjugate 
to?+p Ao (T)-p Ao ($J 1 (x)) along 73,-1 

(S),-iy ,t (£) we ma y m fact obtain s(x, to; 0, t) for 
all t € [?+ Pa (2) — Pa *o] for which it is defined by continuity. Finally since 

7*-MS) -vot o C5r)(( >e + PA (2) - PAoC^J 1 ^)))) C M\M we can also find s(x,t ;0,i) 
for all t S (0, to] for which it is defined. 

Now suppose that t' Q 6 (— e, pa ( < 1 > J 1 (s)))- Provided that t' and ? are not con- 
jugate along 7a; ,r;o: there exists A € 7^ j£ such that z S domain(p A j ) ). We may 
now repeat the previous construction with to replaced by t' to obtain s(x, t' o ;0,t) 
for all t e (0, t ] where it is defined, and this gives all of the data necessary to ap- 
ply Theorem |1.1| Therefore we can recover the wave speed v in a neighborhood of 
-v x (x) ((^' *o + e ))- By the hypotheses W can be covered by sets of this form, 
and so we can recover the wave speed v on all of W' as claimed. □ 

7. Conclusion. We generalized the method of Dix for reconstructing a depth 
varying velocity in a half space, where depth is the Cartesian coordinate normal to 
the boundary, to a procedure for reconstructing a conformally Euclidean metric on a 
region of R™ from expansions of diffraction travel times generated by scatterers in the 
region and measured on its boundary. Our procedure consists of two steps: In the first 
step, we reconstruct the directional curvature operator along geodesies as well as the 
metric in Riemannian normal coordinates. Riemannian normal coordinates can be 
thought of as "time" coordinates as they appear in so-called seismic time migration. 
We note that the directional curvature operator did not appear in the method of 
Dix because of the class of velocity models he considered. In the second step, the 
velocity and the geodesies on which the velocity is reconstructed are obtained through 
a transformation to Cartesian coordinates; this can be thought of as a generalization 
of the "time-to-depth" conversion in the framework of Dix' original formulation. In 
dimension three or more both steps are essentially formulated in terms of solving a 
closed system of nonlinear ordinary differential equations, for example, by application 
of a Runge-Kutta method. In dimension two the second step requires the solution 
of a Cauchy problem for an elliptic operator which may suffer from stability issues. 
Through the associated discretization, we accommodate the case of a finite number 
of scatterers in the manifold. We admit the formation of caustics. 

Appendix A. Converting travel times measured at the boundary to the 
shape operator. 

In this appendix we show more explicitly how the somewhat abstract procedure 
described in the proof of Theorem |1.2| to move from travel times, or generalized 
distance functions, to the shape operators for wavefronts can be done in a particular 
case. We assume that M is a lower half space and v is constant in a neighborhood of 
the boundary and thus can be extended smoothly as a constant. Here we have the 
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Fig. A.l. On the left is illustrated a few rays originating at a common diffraction point and 
contributing to a single wavefront. The travel times are measured at z = producing a curve in 
(x,t) space as illustrated on the right. 



application in mind, and so we use the term travel time for the distance in the metric 
9- 

We use Cartesian coordinates (x,z) € W 1 ^ 1 x R and assume that we have a flat 
boundary at dM = {z = 0} so that M = {z < 0} and M = R™. We focus on a 
single diffraction point, yo> arL d measure the travel times from yo as a function of the 
transverse coordinates at {z = 0} which we label as x. We assume that the wave 
speed v is equal to a constant c in a neighborhood of {z = 0} and write the travel 
time as t which is a function of x for x in some open set U C 

First we can determine the tangent vector of the ray passing through a certain 
point x of the boundary and corresponding to the travel time t(x) by 

dz ^ / 

: = -^:(t(x),x) = \fc 2 - \x(x)\ 

(these equations should be compared with ( |6.2[ )). We are using the prime notation here 
to indicate that the indices only run from 1 to n — 1. In practice it would be possible 
to calculate x by finding the normal to the surface (x,t(x)) which is proportional to 
the vector 

(-x(x),c 2 ) . 

Also note that if we look at multiple wavefronts coming from different diffraction 
points, then we can identify contributions from diffraction points along the same ray 
by the fact that they must satisfy 

±i(x) = x 2 {x) 



(compare with (6.1)) where i\ and xi correspond to the two different wavefronts. 

Now we choose to > sup t and consider the surface of points with travel time to 
from the given diffraction point. In the body of the paper this surface is labeled So,t - 
We may parametrize the surface So,t by x using the map (labeled in the body 
of the paper) 

x t-^ (x = (to — t(x))x(x) + x, z — (to — t(x))z(x)j . 
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It is through this map that we define the coordinates on £ ,t • F° r < r < t — t(x) 
we can also write a formula for the inverse of the coordinate map for the coordinates 
(x, r) (\I> to in the body of the paper). This is 

(x, r) h-» (x = (to — t(x) — r)x(x) + x,z = (to — t(x) — r)z(x)) . 

Now we may calculate the matrices F k and Jj° which give respectively the parallel 
translated fields Fj and the coordinate vector fields d/dx^ (which are also Jacobi 
fields) in the Cartesian coordinates (note that J k differs from jj° in the body of the 
paper which represents the Jacobi fields with respect to the parallel frame). To take 
advantage of the summation notation we also write x n = r and x n — z. Then we 
have 



d 



dx k d 



(x,z) 



where we may calculate (using the notation s(x,r) = t — t(x) — r) 



3 c 2 

—x k ' —z 



The parallel fields are the same as the coordinate vectors at r = 0, and since for 
< r < to — t(x) the wave speed is constant they are in fact constant in the Cartesian 
frame on the same set. Thus we have, for < r < t — t(x), 



F(x,r,t ) = I 1 y k ' Z 



These matrices may also be written directly in terms of the travel time function t. 
Indeed if we write Hess(i) for the Hessian matrix of t, and Vi for the gradient then 

(s + sc 2 Hcss(t) - c 2 (V*)(V*) T - , sc2 - , Hess(t)Vt- c 2 J c~ 2 - |Vt| 2 V 
J(x,r,t ) = V^ 2 -I v *l 



and 



(s + (t - t)c 2 Hc !a (t) - c 2 (Vt)(Vt) T - / t 0- f > c H.B,(i)Vt - c 2 ,/c-2 - |Vt| 2 V 
P(2,r,t )= ^ja- 2 -\Vt\ 2 

\ -a 2 (Vi)T -c 2 ^/c~2 _ | V? |2 

Now we calculate the matrix s(x, t ; , r, to) which corresponds with the shape operator 
of XV ito with respect to the parallel frame. To begin, recall that 

s k (x,t ;,r,t ) = (/*(£, O.S^F^.r)) 

where {/ fc }J? =1 is the dual frame for {fj}™ =1 . Now we can use the fact that the 
Christoffel symbols in Cartesian coordinates are zero where the wave speed is constant 



Reconstruction of a Riemannian metric 



21 



to calculate 

pi I dx^_d i_ dz d\ 

3 \ dxt dxP' dxi g z J 



Therefore 




Once again, this can also be changed to an expression in terms of t and its derivatives 
as 

„ , /c 2 Hess(t) . c2 , Hess(t)Vt\ 

s k 3 {x 1 t,-rM=FJ- 1 \ q W v^^f jF- 1 

and we may use this formula to calculate Sj(x, to; to — t(x),to) which we invert to give 
the initial conditions for the system ( |3.1[ ). To calculate the derivatives with respect 
to t we must consider multiple wavefronts, and using the comment above identify 
singularities in the wavefront corresponding to diffraction points along the same ray. 
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